Ordering phenomena in cooling granular mixtures 
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We report two phenomena, induced by dynamical correlations, that occur during the free coohng 
of a two-dimensional mixture of inelastic hard disks. First, we show that, due to the onset of velocity 
correlations, the ratio of the kinetic energies associated with the two species changes from the value 
corresponding to the homogeneous cooling state to a value approximately given by the mass ratio 
m\/m2 of the two species. Second, we report a novel segregation effect that occurs in the late stage 
of cooling, where interconnected domains appear. Spectral analysis of the composition field reveals 
the emergence of a growing characteristic length. 
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A gas of inelastic hard spheres (IHS), due to its rela- 
tive simplicity, represents a standard reference model for 
fluidized granular materials 0. Many of its properties 
are well understood. In particular, a great deal of atten- 
tion has been devoted to the study of the cooling process 
that occurs when an assembly of grains, initially in mo- 
tion, evolves in the absence of any external energy feed. 
In this case, while the asymptotic state (all particles at 
rest) is trivial, the dynamics of the cooling process dis- 
plays several interesting features. Three stages can be 
identified: 

i) After an initial transient, during which the velocity 
distribution is not necessarily Maxwellian , the system 
enters the homogeneous cooling state (HCS), where the 
average kinetic energy decreases according to Haff 's law 
i, K{t) = K{Q)l{l+tltof. Such a state is characterized 
by a uniform density and a lack of strong correlations 
among the velocities 0. 

ii) An intermediate regime, wherein the system re- 
mains homogeneous, but the velocity field develops vor- 
tices (fig. These vortices result from the paralleliza- 
tion of velocities induced by inelasticity. Due to the pro- 
gressive build up of velocity correlations, the kinetic en- 
ergy decays as t"''/^ |3,IE0- 

iii) Goldhirsch and Zanetti suggested that vortices 
represent the physical mechanism leading to the third 
stage, wherein the density field itself loses its homogene- 
ity and dense clusters appear, surrounded by rarefied re- 
gions. At this stage neither the velocity nor the temper- 
ature fields are homogeneous anymore (see fig.Ql. 

Most of the existing information on the cooling proper- 
ties of granular matter concern one-component systems. 
A few investigations on binary mixtures |^, focused 
on the HCS only. Garzo and Dufty |3| showed that for 
a binary mixture the average kinetic energies per par- 
ticle - i.e. the kinetic temperatures Ki and K2 of the 
two species - are different, and their ratio 7 = K1/K2 
remains constant in the HCS. 

In this Letter we explore the cooling dynamics of a mix- 
ture of inelastic hard disks across the previously listed 
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FIG. 1: Status of the system at points A (upper panel) and 
B (lower panel) of fig. |21 The instantaneous velocity field is 
superimposed. Early vortices can be seen in the upper panel, 
with particles still (approximately) uniformly distributed. In 
the lower panel, dense clustering has appeared and the typical 
size of vortices has significantly increased. Strong correlations 
in the velocity field are visible. 



dynamical regimes by using an event-driven molecular 
dynamics simulation and a recently introduced lattice 
model [llj . We show that during the second stage of cool- 
ing - at the onset of local velocity correlations - the ra- 
tio of the kinetic temperatures of the two species sharply 
departs from the nearly steady value attained during the 
homogeneous cooling stage, and eventually approaches 
an asymptotic value that only depends on the mass ratio 
of the two species. In the last stage we find evidence for 
segregation of the species, which is revealed both by lo- 
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FIG. 2: Ratio of the kinetic energies of the two species for 
different values of the coefficient of restitution. The crossover 
from K1/K2 ~ 1 to K\/K2 ~ m\/m2 is clearly visible. For 
comparison, we also show the corresponding result for the 
lattice model. The status of the system at points A and B is 
shown in fig.0 

cal indicators and by spectral analysis of the composition 
field. 

For our model we consider a two-dimensional (2D) sys- 
tem containing A^i = 5 x 10^ and A''2 = A''i disks of equal 
radius R and masses mi and 7712, respectively. Initially, 
the discs are uniformly distributed in a square domain 
of side L — 10"^ i?, with periodic boundary conditions, 
and their velocities are set by sampling Maxwell distri- 
butions corresponding to the same temperature. In the 
IHS model the collisions are instantaneous, conserve total 
momentum and dissipate a fraction of the kinetic energy. 
The velocity vj of a particle j after collision is 

Vj = Vj + /iy [1 + r(w„)] w„ n , (1) 

where n = (r^ — )/ jr^ — | is the unit vector joining the 
centers ri and rj of the disks at contact, u„ = (vi — vj) -n 
is the relative normal velocity of the disks before collision, 
/i.y = mi/{mi + rrij), and r(w„) is a velocity-dependent 
coefficient of restitution. Inelastic collapse is avoided by 
choosing a r(vn) such that, for relative normal velocities 
falling below a threshold vq, r approaches unity when 
Vn — + 0, according to the power law given in ref. 0. 
Different cutoff velocities were simulated in order to as- 
certain the robustness of the investigated properties. An 
efhcient event-driven simulation allowed us to probe the 
late stages of the coohng process ■ 

During the initial homogeneous cooling phase, the be- 



havior of the mixture is well described by the Boltzmann- 
Enskog (BE) transport equation |9(|. Theoretical values 
for the ratio 7 can be computed within this framework, 
and their values for the chosen set of microscopic param- 
eters are shown in fig. |21 (dashed lines). These theoretical 
values agree with the initial plateau observed in the tem- 
poral evolution of 7. However since the BE approach, 
by construction, neglects velocity correlations, one can 
expect that the intermediate cooling regime cannot be 
properly described within its framework. Indeed, hydro- 
dynamics predicts that the HCS is unstable with respect 
to long-wavelength fluctuations of the velocity field. Its 
repercussions on the behavior of a mixture have not been 
explored to date. We verified that Ki{t) and K2{t) ex- 
hibited - separately - the t^^ and t"^ scaling behaviors 
that are commonly observed in one-component systems. 

Fig. El shows that 7 = K1/K2 changes rather sharply 
from the initial plateau value of the HCS to a higher 
plateau value. The number of collisions at which the 
crossover occurs is proportional to (1 — r'^)~^ and di- 
verges when the system approaches perfect elasticity. As 
we mentioned above, the position of the lower plateau 
is well accounted for by the BE theory. From fig. |21one 
can also see that the asymptotic, higher plateau corre- 
sponds to K1/K2 — mi/m2. This correspondence was 
confirmed for several values of the density ratio of the 
two species. Whereas in the HCS the inelasticity deter- 
mines only a small departure from equipartition, in the 
intermediate and late regime such a breakdown becomes 
more severe: the specific asymptotic value of 7 means 
that particles with different masses eventually share the 
same velocity distributions, i.e. "equipartition" of veloci- 
ties eventually sets in. Although not previously reported, 
the phenomenon has a simple explanation: due to the in- 
elasticity of collisions, the velocities of the particles tend 
to become more and more parallel. Since momentum con- 
servation does not allow the velocities to select a unique 
direction, the velocity field develops many vortices (see 
fig.^, i.e. a situation of local alignment. The occurrence 
of a collision depends only on the relative velocities of 
particles, and not on their masses, coefficients of resti- 
tution or other properties. Therefore, the velocity of a 
given particle tends to be locked to the velocity of neigh- 
bors, regardless of their species. The fact that the ratio 
K1/K2 eventually attains the value mi/m2 is a striking 
manifestation of the long range velocity correlations that 
develop in a granular fluid. We checked the above argu- 
ment against a different case, wherein densities are kept 
unchanged and the masses of the two species are made 
equal by choosing different radii for them. The result is 
shown in fig. (gray diamonds). In this case the value 
of 7 stays close to mi/m2 = 1. Also, the dependence 
on other parameters (such as different restitution coeffi- 
cients, geometrical factors, tangential friction modeled as 
in ref. [IjI) was found to be rather weak. Of course, real 
granular materials are subject to more complex dissipa- 
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tive forces. Our tests, however, show that 7 asymptot- 
ically depends mainly on mi/m2, although the detailed 
crossover dynamics reflects the specific parameters. 

It is interesting to observe that the previously pre- 
sented scenario is reproduced by a very simple extension 
to mixtures of the inelastic lattice model of ref. 
The mixture consists of A^i and N2 "particles" with the 
same physical properties as above, sitting randomly on 
the M = iVi + A^2 nodes of a 2D triangular lattice. Every 
particle is characterized by a 2D vector velocity Vi, but its 
position stays fixed. A dynamics is induced by randomly 
choosing two neighboring particles and updating their 
velocities according to the rule of eq. The free stream- 
ing component of the dynamics is neglected, but velocity 
correlations are properly taken into account. After an 
exponential decay, the energy decreases as if ~ r"''/^, 
as usual. Correlations become more and more important 
as time elapses, showing the presence of a growing length 
scale L(t) ~ r-'^/^ associated with the average diameter 
of the vortices As shown in fig. |21 (gray triangles) 

the lattice mixture model can reproduce the asymptotic 
value of 7. To demonstrate that velocity correlations are 
responsible for the specific value of 7, we shuffled the ve- 
locity of each particle right after collision, by exchanging 
its velocity with that of a randomly chosen particle. This 
suppresses completely the t~'^I'^ energy decay, the vortex 
instability, and the crossover in 7. 

In the third stage of cooling, density fluctuations 
lead to cluster formation, like in mono-disperse sys- 
tems. In addition to that (see fig. 0)), we observe a 
novel and appreciable segregation of species, that we 
quantify by studying the following ratios S\ and 6*2 of 
integrated (partial) pair distribution functions: Sa = 
/o 9aa{r)rdr/ gap{r)rdr, where ga(i(r) {a, (3 = 
1, 2) are the partial pair correlation functions and Z) is a 
cutoff distance, here chosen equal to 5R. Fig.^shows the 
time evolution of and ^2 for different choices of the 
mass ratio. For values of r corresponding to the transi- 
tion in 7 (refer to fig. El, local segregation appears: heavy 
particles tend to group together (5*1 > 1) and the effect 
increases with the mass ratio. A simple heuristic argu- 
ment can explain such a finding. The average change of 
kinetic energy {AEa}p of a particle of species a colliding 
with a particle of species (3 can be expressed as a function 
of the local kinetic temperatures Ti and T2 as [l3^ 
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One sees that particles with a smaller mass gain energy 
((A£^2)i > 0) at the expenses of the heavy particles, 
provided that T1/T2 > ((1 - r)Tni/m2 + 2)/(l + r) > 1. 
Indeed, our numerics confirm that this condition is satis- 
fied for our system. For a perfectly elastic system there 
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FIG. 3: Local segregation; the temporal evolution of the 
ratios 5*1 and S2 of integrated correlation functions (see text) 
is shown for m\/m2 = 2 and m\/m2 — 8. To assess noise 
levels, a control run with mi = m2 is also shown. 



is no energy transfer, since Ti = T2. In the dissipative 
case, instead, since the ratio of the kinetic temperatures 
T1/T2 > 1, light particles - on average ~ gain kinetic 
energy when colliding with heavy particles. This re- 
sults in an effective repulsion of the light particles by the 
heavy ones. The above scenario is also consistent with 
the asymmetry between species observed in the temporal 
evolution of the indicators Sg^p of local segregation. 

Such a microscopic mechanism driving segregation has 
also global effects. Let us consider the long- wavelength 
properties of the composition field (/)(r,i) = ni(r,i) — 
n2{r,t), the difference of the local number densities of 
the two species, whose dynamics is slow due to its con- 
servative nature. The analysis of the structure factor 
F(k,t) = 4)(k,t)4)*(k,t) (where (^(k,t) is the Fourier 
transform of (/)(r,t)) reveals a remarkable scenario: at 
the onset of clustering the uniform composition state be- 
comes unstable and evolves into an interconnected pat- 
tern of domains with different composition, as shown in 
fig. 21 Long- wavelength fluctuations of (f>{r, t) are am- 
plified and - contrary to the usual diffusion processes - 
particles of different species tend to segregate. The tem- 
poral evolution of the above (circularly averaged) struc- 
ture factor is reported in fig. [S] In the homogeneous 
stage of cooling, the structure factor is featureless (right 
triangles). When clustering appears, a broad peak de- 
velops (diamonds) and with time it grows higher, while 
shifting towards low wavevectors. Reshuffling the parti- 
cle species - while keeping their positions unchanged - 
suppresses the peak altogether, as shown in the inset of 
fig- El This proves that the above spectral properties are 
truly a manifestation of progressive, long-wavelength or- 
dering of the two species. Moreover, the peak in F{k, t) is 
unrelated to the short wavevector peak that - because of 
clustering - is expected in the structure factor of the to- 
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FIG. 4: Composition pattern within a cluster after r — 1300 
collisions per particle, for r = 0.95 and m\/m2 = 8. Black 
discs correspond to heavier particles. The radius of the par- 
ticles is R = 0.1 cm. In inset: global state of the system. 

tal density field (see fig. El inset): the composition field 
and the density field have very different characteristic 
lengths. 

In conclusion, we found that the strong velocity cor- 
relations induced by inelasticity result in a ratio of ki- 
netic energies K1/K2 — mi/m2 that is not accounted 
for by the BE theory, which underestimates the effect of 
dynamical correlations. This can be regarded as a man- 
ifestation of the noise reduction effect 0| of granular 
systems, where a suitable collective arrangement of the 
velocities minimizes dissipation. In addition to that, in 
the clustered stage, segregation builds up according to 
a phenomenology which bears similarities to the spin- 
odal decomposition occurring in metallic alloys, where a 
growing characteristic length appears in the composition 
field. C. C. gratefully acknowledges support from the 
Japan Society for the Promotion of Science (JSPS). 
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FIG. 5: Circularly averaged composition structure factor, 
F(k,t), at different times. With time, the peak grows higher 
and moves towards smaller wavevectors. The inset shows how 
the peak for r = 1300 (solid circles, right vertical axis) is sup- 
pressed on reshuffling the particle species (empty circles, right 
axis). The structure factor of the total density field (empty 
squares, left axis) has a peak at a much lower wavevector. For 
these curves, r = 0.95 and m\/m2 = 8, as in fig.|l] 




